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Coalescing binary black hole mergers are expected to be the strongest gravitational wave sources 
for ground-based interferometers, such as the LIGO, VIRGO, and GE0600, as well as the space- 
based interferometer LISA. Until recently it has been impossible to reliably derive the predictions 
of General Relativity for the final merger stage, which takes place in the strong-field regime. Recent 
progress in numerical relativity simulations is, however, revolutionizing our understanding of these 
systems. We examine here the specific case of merging equal-mass Schwarzschild black holes in 
detail, presenting new simulations in which the black holes start in the late inspiral stage on orbits 
with very low eccentricity and evolve for ~ 1200 M through ~ 7 orbits before merging. We study 
the accuracy and consistency of our simulations and the resulting gravitational waveforms, which 
encompass ~ 14 cycles before merger, and highlight the importance of using frequency (rather 
than time) to set the physical reference when comparing models. Matching our results to PN 
calculations for the earlier parts of the inspiral provides a combined waveform with less than half a 
cycle of accumulated phase error through the entire coalescence. Using this waveform, we calculate 
signal-to-noise ratios (SNRs) for iLIGO, adLIGO, and LISA, highlighting the contributions from the 
late-inspiral and merger-ringdown parts of the waveform which can now be simulated numerically. 

Contour plots of SNR as a function of z and M show that adLIGO can achieve SNR > 10 for some 
IMBBHs out to z ~ 1 , and that LISA can see MBBHs in the range 3 x 10 4 < M/M © < 10 7 at 
SNR >100 out to the earliest epochs of structure formation at z > 15. 

PACS numbers: 04.25. Dm, 04.30.Db, 04.70. Bw, 04.80.Nn 95.30.Sf, 95.55.Ym 97.60.Lf 


I. INTRODUCTION 

Coalescing binary black holes (BBHs) are strong 
sources of gravitational radiation for ground-based de- 
tectors such as LIGO, VIRGO, and GEO600, and for the 
space-based LISA. This coalescence is generally consid- 
ered to proceed in three phases: inspiral, merger, and 
ringdown [1]. During the inspiral stage, the black holes 
are well separated and spiral together on quasicircular or- 
bits. This is followed by the merger stage, during which 
the black holes plunge towards the center and merge to- 
gether, forming a common horizon. This distorted rem- 
nant then ‘rings down’ to form a Kerr black hole by shed- 
ding its nonaxisymmetric modes as gravitational radia- 
tion. 

Different techniques are used to calculate the gravi- 
tational radiation from these stages. The inspiral can 
be treated analytically with post- Newtonian (PN) tech- 
niques; the gravitational waveform is a chirp, a sinusoid 
that increases in both amplitude and frequency. The 
later part of the ringdown can be handled analytically 
using black hole perturbation theory, with the result- 
ing waveforms being exponentially damped sinusoids of 
constant frequency. However, calculating the waveforms 
from the dynamical merger, which produces the highest 


luminosity signal, requires numerical relativity simula- 
tions of the full Einstein equations in three spatial di- 
mensions plus time. With the first generation of ground- 
based detectors now taking data and LISA moving for- 
ward, knowledge of the merger waveforms and their im- 
pact on detectability and parameter estimation is urgent. 

Recently there has been dramatic progress in numeri- 
cal relativity calculations of BBH mergers. The first full 
orbit of an equal mass nonspinning BBH was achieved 
nearly three years ago [2] (see also [3]), using a confor- 
mal formalism and comoving coordinates, with the black 
holes represented as punctures [4], and held fixed in the 
grid. This was followed about a year and a half later 
by the first simulation of the final orbit, merger, and 
ringdown [5]; this work was carried out using general- 
ized harmonic coordinates with excised black holes mov- 
ing through the computational domain [6, 7]. Roughly 
six months later the moving puncture method, which is 
based on a conformal formalism and allows the puncture 
black holes to move across the grid without the need for 
excision [8, 9], was introduced. This simple yet power- 
ful technique proved to be highly effective [10, 11, 12], 
allowing the evolution of black holes through multiple 
orbits, merger, and ringdown [13, 14]. Simulations with 
unequal masses [15, 16, 17] and with spins [18, 19] quickly 
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followed. Longer simulations with several orbits before 
merger have also been carried out using excision and har- 
monic techniques [20] , allowing for comparisons between 
the resulting waveforms [21]. 

This explosion of new work on BBH mergers has 
opened up exciting opportunities for applications of these 
results. A number of key applications center on gravita- 
tional wave data analysis. In particular, the need for 
accurate templates for the full gravitational wave sig- 
nal - encompassing the inspiral, merger, and ringdown 
phases - means it is essential to understand the relation- 
ship between PN waveforms and the results from numer- 
ical relativity [14, 20, 22]. And, with longer simulations 
now becoming available, it is becoming feasible to apply 
these waveforms to questions about the detectability of 
BBHs for various detectors [20]. In this paper, we focus 
especially on issues of accuracy in long simulations, start- 
ing in the late-inspiral regime, and their applications in 
gravitational wave data analysis. 

The most rapid advances cover the previously least un- 
derstood portion of the radiation, the final few cycles 
of radiation generated from near the “innermost stable 
circular orbit” (ISCO) and afterward which we call the 
“merger-ringdown” . Already there has been considerable 
progress toward a full understanding of this important 
portion of the waveform through which the frequency 
sweeps by a factor of. ~ 3 up to ringdown. A signifi- 
cant development in this regard was the demonstration of 
initial-data-independence of merger waveforms for equal- 
mass, spinless black holes [14]. Today, there is general 
agreement that the merger of such a system produces a 
final Kerr BH remnant with spin a/m ~ 0.7, and that the 
amount of energy radiated in the form of gravitational 
waves, starting with the final few orbits and proceeding 
through the plunge, merger and ringdown, is E t&( i ~ 4% 
[5, 8, 9, 13, 14]; see [23] for a review. There is also 
consensus on the overall simple shape of the waveforms; 
detailed comparison of results among numerical relativ- 
ity groups has already begun [21], with more to follow. 
The exploration of parameter space, including various 
mass ratios [15, 16, 17] and spins [18, 19], is now un- 
derway. Future efforts will involve pushing this frontier 
to increasingly complex mass-ratio and spin-orientation 
combinations, and to establishing initial-data indepen- 
dence across these parameters. 

Though more challenging, it is now becoming practical 
to simulate BBHs starting in the late inspiral as well. We 
report on new simulations covering roughly an additional 
factor of 3 in frequency before the “merger-ringdown”. 
We have simulated the coalescence of equal-mass, non- 
spinning black holes, starting in the late-inspiral regime, 
and continuing through the merger and ringdown. The 
black holes start on orbits with very low eccentricity and 
spiral together through ~ 7 orbits before merging. The 
resulting gravitational waveform has ~ 14 cycles, en- 
abling a robust consistency test with late inspiral PN 
waveforms as reported in [22]. In Sec. II we present an 
overview of our numerical results. Readers with an in- 


terest in the details of the simulations will find them in 
Sec. Ill, where we discuss our methodology and the accu- 
racy of our simulations, as well as introduce the impor- 
tant notion of considering physical quantities as functions 
of frequency rather than time when comparing predic- 
tions from different computations. 

The rapid advances in understanding merger-ringdown 
waveforms with emerging simulation results for the late 
inspiral create a new context for considering gravitational 
wave observations of this part of the signal from such a 
system. In Sec. IV we revisit an analysis of merger de- 
tectability originally carried out in Flanagan and Hughes 
[1], now in the context of present-day and emerging 
waveform modeling capabilities. In particular, we ex- 
amine the relative contributions to the signal-to-noise 
ratio (SNR) for the merger-ringdown and late-inspiral 
portions of equal-mass coalescence signals, as observed 
by initial LIGO (iLIGO), advanced LIGO (adLIGO), or 
LISA. We also plot contours of SNR as functions of red- 
shift z and total mass for adLIGO and LISA, showing 
their ability to detect BBHs in the cosmos. These SNR 
calculations draw on a refined waveform model based on 
a best-estimate waveform utilizing both numerical rela- 
tivity and PN results. We conclude with Sec. V, which 
contains a summary and discussion of our main results. 


II. OVERVIEW OF SIMULATION RESULTS 

Late-inspiral evolutions, covering more than a few or- 
bits prior to ringdown, are more challenging than shorter 
merger-ringdown simulations. In this regime there is a 
stronger requirement for numerical stability and a greater 
demand for computational resources. In addition, better 
accuracy is required to control the accumulated phase 
error, which in turn constrains the numerical error that 
can be tolerated in the rate of energy loss. 

Using the moving puncture technique [8, 9, 10] we 
have simulated the evolution of a nonspinning equal-mass 
BBH starting at relatively wide separation, ~ 1200M or 
~ 7 orbits before the formation of a common event hori- 
zon. Here M is the total mass the system would have 
had when the black holes were very far apart and be- 
fore radiative losses became significant. We used fourth- 
order-accurate finite differencing techniques with adap- 
tive mesh refinement (AMR) to resolve the dynamics near 
the black holes and in the wave propagation region. We 
carried out three runs using similar grid refinement struc- 
tures, but at different resolutions: low ( hf = 3M/64), 
medium (hf = 3M/80) and high (hf = M/32). Here, hf 
is the grid spacing in the regions with the highest reso- 
lution in each simulation, those being the regions around 
each black hole. Overall, the Hamiltonian constraint con- 
verges at fourth order, and the momentum constraint at 
better than second order, throughout the runs. 

Fig. 1 shows the resulting gravitational waveform in 
the context of recent developments in black hole evolu- 
tions. The solid (blue) curve gives the gravitational wave 
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FIG. 1: Gravitational strain waveforms from the merger of equal-mass Schwarzschild black holes. The late part of the merger 
(t > —50 M) is robustly determined and relatively easily calculable, while simulations of the late inspiral (early part of the 
waveform) are rapidly approaching the phasing accuracy required for observational applications [22] . The solid blue line shows 
our current “best” numerical waveform. The dashed red line shows a comparison waveform from a run starting with the same 
initial data as R4 in Ref. [14] and the dash-dotted green curve shows the results from the highest resolution R1 run in Ref. [14]. 
All waveforms have been extracted at R e xt = 40M and shifted in time so that the moment of maximum 1/4 radiation amplitude 
occurs at time t = 0. 


strain from the dominant l = 2, m = 2 spin-weighted 
spherical harmonic from our highest resolution simula- 
tion, extracted at R ext = 40M. This represents an ob- 
servation made on the equatorial plane of the system, 
where only a single polarization component contributes 
to the measured strain. Time t — 0 is taken to be the 
moment of peak radiation amplitude as measured by the 
Weyl curvature ^4; the symbols along the time axis mark 
the points at which the system reaches the ISCO calcu- 
lated for a Schwarzschild black hole (circle), EOB 3PN 
[24, 25] (diamond), and adiabatic 3 PN [25, 26] (square). 
For comparison, we show two other waveforms from ear- 
lier simulations carried out by this group; both were ex- 
tracted at jR ext = 40M and have been shifted in time 
and initial phase so (as in [14]) that the moment of peak 
t />4 radiation amplitude occurs at t = 0. As these dif- 
ferent simulations may radiate different fractions of the 
initial mass, we choose the mass m/ of the post-merger 
Kerr hole as the natural length scale for comparison (see 
discussion of Table I for details). Because of radiative 
losses, rrif « 0.95 M. 

The dashed (red) curve shows the results from a model 
that starts ~ 550 M before merger with the same initial 
data as run R4 in Ref. [14] but using a gauge more 
conducive to numerical accuracy [10]). The dot-dash 
(green) curve shows a simulation that starts ~ 200 M 
before merger; this is the high resolution run R1 from 
Ref. [14]. All three waveforms agree to within ~ 1% for 
the merger-ringdown burst part of the waveform, starting 
near t ~ —50 M. 

Astrophysically, BBHs in this relatively late stage of 
their evolution are expected to be traveling on nearly 
circular orbits, as any initial eccentricity would have 


been radiated away early in their evolution. In this new 
model, we start the black holes on nearly circular orbits 
with very small eccentricity e < 0.01, as defined below. 
The resulting black hole trajectories are shown in Fig. 2, 
where the tracks mark the paths of the punctures; for 
clarity, only a single black hole from each simulation is 
shown. The dashed (blue) curve shows the results from 
our new high resolution long run; note that the early 
orbits are very nearly circular. The solid (red) curve 
shows the trajectory of the moderately long comparison 
run (shown by the dashed line in Fig. 1). At early times, 
this model is significantly more eccentric than our new 
results. At later times this trajectory locks on to that of 
our new run, ~ 2.5 orbits before merger. 

The black hole separation as a function of time is 
shown for these two models in Fig. 3. The greater ec- 
centricity of the previous model (solid curve) is clearly 
distinguished here. We also show all three resolutions of 
our newest model. The slight deviations from the over- 
all smooth trend give an indication of the small amount 
of eccentricity in these latter simulations. Note however 
the differences between these three resolutions, particu- 
larly in the total time between the start of the simulation 
and the merger. Although significant, the differences in 
merger time appear to converge at a rate consistent with 
fourth order error. 

In the next section we consider numerical techniques 
and the fidelity of the simulations in more detail, includ- 
ing the convergence and conservation properties. We also 
carry out a detailed analysis of the resulting gravitational 
waveforms, underscoring the differences between using 
time and frequency to set references for comparing mod- 
els. Readers more interested in issues of detectability and 
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FIG. 2: The trajectory of one of the binary system’s black 
holes through ~ 7 revolutions before coalescence for our high 
resolution case is shown by the dotted line. The solid line 
gives the trajectory of the moderately long comparison run. 



FIG. 3: The coordinate separation between the puncture 
black holes is shown as a function of time. The solid line shows 
the results for the comparison run, which has relatively large 
eccentricity. The other curves show the three resolutions for 
our new simulations, all having noticeably less eccentricity. 


SNR analyses for iLIGO, adLIGO, and LISA are invited 
to skip to Sec. IV. 


III. SIMULATION DETAILS 

A. Numerical Methodology 

For initial data we used Brandt-Briigmann puncture 
data[4], generated by a second-order- accurate multigrid 


solver, AMRMG[27]. The puncture parameters were deter- 
mined by the Tichy-Brugmann prescription for quasi- 
circular initial data[28], adjusted slightly, as informed 
by previous empirical experience, to reduce eccentricity. 
The success of this approach in giving a circular inspi- 
ral is roughly indicated by the puncture track (dotted 
line) in Fig. 2 and the curves in Fig. 3. The punc- 
ture track was computed by integrating the equation 
Xpunc = —0(x pU nc), where x punc is the position of the 
puncture, and (3 the shift. 

This data was evolved using standard BSSN evolu- 
tion equations, with the addition of dissipation terms as 
in [29] and constraint-damping terms as in [30] in or- 
der to ensure robust stability. The gauge condition was 
that recommended in [10] for moving punctures. Time 
integration was performed with a fourth-order Runge- 
Kutta algorithm, and spatial derivatives with fourth- 
order-accurate finite differencing stencils. For the outer 
boundary we employed a second-order-accurate Sommer- 
feld condition, pushed to \x{\ = 1536M to keep reflections 
far from the source. Adaptive mesh refinement was im- 
plemented via the software package PARAMESH [31 , 32] , 
and interpolation between refinement regions was fifth- 
order accurate. 


B. Simulation Analysis 

The accuracy of the simulations was assessed by var- 
ious means. We first considered the L\ norm of the 
constraints in each refinement region, the grid structure 
having been designed to be commensurate for all resolu- 
tions; results for the finest (top panel) and second finest 
(bottom panel) refinement regions are plotted in Fig. 4 
and Fig. 5. The Hamiltonian and momentum constraints 
were found to be convergent at an apparent order of 2.5 
in the finest grid, where the error at the puncture can be 
expected to dominate and fourth-order finite differenc- 
ing must break down due to the irregularity there. In all 
coarser regions the Hamiltonian constraint appears to be 
fourth-order-convergent, while the momentum constraint 
appears closer to second-order convergent. 

From each simulation we have measured the radia- 
tion in the form of the complex Weyl tensor compo- 
nent 0 4, specified consistently with [35] to leading order 
in 1/r. The gravitational wave strain is related to 04 
by /i + + ih x = 2-04, and can be recovered by integra- 
tion, with the two complex integration constants chosen 
to keep the strain close to oscillating about zero. For 
some applications we also examine waveforms in the form 
of the strain rate v = h+ + ih x , which is the quan- 
tity from which radiation energy is directly obtained. 
To extract the waveform information from the simula- 
tion we define a series of coordinate spheres of differ- 
ent radii Re Xt /M € {40,60,100} on which we measure 
spin-weighted spherical harmonic components of 04 via 
a second-order algorithm described in [33, 36]. In this 
paper we focus exclusively on the l = 2,m = 2 compo- 
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FIG. 4: Convergence plot for the Hamiltonian constraint Ch- 
The top panel shows results from the finest grid and has been 
scaled so that for 2.5 order convergence the curves should 
superpose. The bottom panel shows results from the second 
finest grid and has been scaled so that for fourth order conver- 
gence the curves should superpose; the curves indeed appear 
to be fourth order convergent. 


FIG. 6: V >4 waveform calculated at three different extraction 
radii, and scaled by the approximate Schwarzschild areal ra- 
dius. Times have been shifted according to the approximate 
Schwarzschild “tortoise coordinate” [33, 34], as appropriate to 
compare the R ex t = 40 M and R ex t = 60 M waveforms with 
the R ext = 100M waveform. 




FIG. 5: Convergence plot for the Momentum constraint Cm- 
The top panel shows results from the finest grid and has been 
scaled so that for 2.5 order convergence the curves should su- 
perpose. The bottom panel shows results from the second 
finest grid and has been scaled so that for fourth-order con- 
vergence the curves should superpose; the curves appear less 
than fourth-order convergent but better than second-order 
convergent. 


nent of the radiation, which mirrors the / = 2,m = —2 
component because of equatorial symmetry. Other com- 
ponents are considerably smaller, containing ~ 1% as 
much energy [14]. 

Fig. 6 compares waveforms from our high-resolution 
simulation extracted on each of our three extraction 
spheres with the i? ex t = 40 M and R ex t = 60M wave- 
forms shifted in time by the intervening propagation time 


derived based on a Schwarzschild black hole background. 
The generally good agreement for all three waveforms in- 
dicates that potentially worrying subtleties related to the 
relatively close distance of the extraction spheres, such 
as non-linear propagation effects, or tetrad-specification 
sensitivity, do not seriously affect the waveforms. On the 
other hand, some small differences are evident. For the 
early portion of the waveforms, shown in the inset, the 
results from the closest extraction sphere show slight dif- 
ferences in amplitude and phase, suggesting that 1/Rext 
radiation details are not yet insignificant here; this is not 
surprising given that the extraction radius in this case 
is only ~ 1/4 of a wavelength. On the other hand, the 
waveform from the farthest extraction radius shows signs 
of dissipation for the later higher-frequency portion of the 
waveform. This is because the radiation has propagated 
significantly farther, through a lower-resolution region on 
the computational grid, by the time it is measured at 
.Rext = 100M. The resolution in most of the intermedi- 
ate region is h — 2 M, about six points per cycle for the 
ringdown radiation. For the rest of the paper, we primar- 
ily employ the waveforms extracted at the intermediate 
distance R ex t = 60M, which is only weakly affected by 
either of the above short- and long- wavelength effects. 

Different information in the waveforms provides inde- 
pendent ways to deduce the energy and angular momen- 
tum of the final black hole produced by the merger. The 
ringdown dynamics of the black hole after merger con- 
tains direct information about the mass and spin param- 
eter of the merged black hole, so that we can deduce m r 
and ( a/m) r by measuring the frequency and decay rate 
of the ringdown radiation. Alternatively, we can mea- 
sure the radiation energy E ra d and angular momentum 
content J ra d at a given radius. Then, since we know the 
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TABLE I: The parameters of the final black hole (mass and 
spin parameter) formed by merger. Here m r and (a/m)r are 
calculated from the complex ringdown frequency u/ r , while 
m/ and (a/m)/ are calculated from the radiation measured 
at the given radius. 


hj 

Rext 

ringdown 
u r m r 

0 a/m) r 

conservation 
771/ (a/771)/ 

3M/64 

40 M 

0.551 - 0.0832 0.971 

0.71 

0.953 

0.70 

3M/80 

40 M 

0.553 - 0.087z 0.945 

0.68 

0.954 

0.72 


60 M 

0.553 - 0.0882 0.931 

0.66 

0.955 

0.70 


100M 

0.553 - 0.0882 0.931 

0.66 

0.959 

0.71 

M/2,2 

40M 

0.553 - 0.0852 0.953 

0.69 

0.955 

0.73 


60M 

0.553 - 0.0852 0.953 

0.69 

0.955 

0.71 


100M 

0.553 - O.O862 0.945 

0.68 

0.958 

0.71 


initial values, conservation of energy and angular mo- 
mentum imply the values of m/ = Madm |t=o — Erad 
and ( a/m)f = (JadmU=o - Jrad)/m 2 . Results from 
both methods are shown in Table I. By comparing rrif 
with m r and ( a/m) r with (a/m)/, we can verify conser- 
vation of energy and angular momentum in the simula- 
tions. In Table I we see that at the highest resolution, 
energy is conserved to within ~ 1% and angular momen- 
tum to within ~ 6%; and the best conservation is seen at 
Rext — 60 M, where energy is conserved to within ~ 0.2% 
and angular momentum to within ~ 3%. 

To an excellent approximation, the l — 2,m = 2 
harmonic of the radiation is polarized in the form ex- 
pected for radiation generated by circular motion. The 
polar representation 0 4 = A^(t) exp(i0^,(£)) of the com- 
plex waveform is particularly natural for circularly po- 
larized radiation, for which A ^ and cj = dcp^/dt vary 
only slowly. The angular polarization frequency u then 
provides a meaningful instantaneous frequency obtained 
directly from the radiation, corresponding to twice the 
angular frequency of orbital motion when the black holes 
are still separate. Because the radiation is measured in 
the weak-field region of our simulations, where gauge de- 
pendence is minimal, this polarization frequency provides 
gauge- invariant information about the binary dynamics. 

If the orbital motion is eccentric, this will leave an im- 
print in the radiation, causing a slight decrease in the 
polarization frequency of radiation generated near apo- 
center. We can recognize eccentric motion by identifying 
periodic deviations from a smooth monotonic trend in 
the time development of the polarization frequency u(t). 

Specifically, we looked at the polarization frequency 
u>(t) = dcp/dt calculated from the strain rate v(t) = 
V(t) exp We see generally similar results whether 
we use strain, strain rate, or 04 to define the frequency, 
but specifically show the strain rate because it comes out 
smoother than 0 4 with respect to small waveform noise, 
but without noticeable detrending issues as in the strain. 
We fit the time dependence of the frequency curves to a 
quartic trend plus an eccentric modulation of the form 


TABLE II: Values resulting from eccentricity fitting. 


Run 

AM 

fl 0 M tlM 2 

$0 

eo 

3M/64 

5 x 10 -4 

0.017 3 x 10 -6 

1 

0.005 

3M/80 

-si 

X 
1— 1 
0 
1 

0.016 4 x 10 -6 

1 

0.007 

M/32 

1 

O 

r— 1 

X 

OO 

0.017 3 x 10 -6 

1 

0.008 


dw = Asin($o + ^o(£ — £o)+^(£ — £o) 2 ) where the quanti- 
ties A, ^0 and tl are fitting constants and to = 61M 
is a time offset approximately accounting for the proper 
gation time to R ext = 60M. Ignoring the early part of the 
simulation where there are transient initial data effects, 
and the late part where the secular trend is very strong, 
we fit the curves over the time range 250 M < t < 1000M. 
We get similar results for eccentricity whether we apply 
high-frequency filtering to the waveform before fitting, 
though we show only low-pass filtered curves in Fig. 7. 
The results of the fitting are summarized in Table II. 

Variations in the details of the fitting procedure, such 
as using strain instead of strain rate or smoothing high- 
frequency noise from c j(t), give results consistent to 
roughly the number of significant digits shown for A and 
H 0, though O and $0 vary more significantly in some 
cases. The period of eccentric oscillations indicated by 
f2o is about 1.5 times the initial orbital period, decreas- 
ing slightly at a rate comparable to the rate at which the 
wave period grows. Allowing A to evolve in time did not 
result in clearly improved fits. 

From these fits we define eccentricity based on the ef- 
fect, in Kepler ian dynamics, of small eccentricity on or- 
bital frequency. Kepler’s second law (conservation of an- 
gular momentum) implies that L oc r 2 u) is constant, lead- 
ing us to define eccentricity as e = A/2u. The constancy 
of A in our fit is consistent with a linear decrease in ec- 
centricity as frequency increases. The initial eccentricity 
eo, obtained in this way, is also given in Table II. 

Even more interesting than the estimates for eccen- 
tricity provided this way are the residual curves u c (t) = 
u(t) — du(t) of angular frequency vs. time with the ec- 
centric part subtracted out. Note that the strictly si- 
nusoidal nature of our fit to the eccentric modulations 
represented by du(t) implies that the underlying secular 
trend should be preserved. Fig. 7 shows the results for 
each of our three simulations together with the unmod- 
ified frequency u(t) for the high-resolution case. The 
plotted curves have also been filtered to remove some 
high-frequency simulation noise present in the early part 
of the simulations. The smooth, now monotonic, trend 
in the curves for u c (t) is an indication that most of the 
wiggles evident in uj(t) were consistent with our model 
for eccentric modulations du(t) ) though the early part is 
affected by transient features related to the shortcomings 
of the initial data model. 

We can take u> c (t) as providing a record of the “harden- 
ing” process, as radiative losses bring the system through 
tightening orbits. The key effect of numerical simulation 
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FIG. 7: Angular frequencies with eccentricity removed, 

aligned such that the frequencies agree when the hf = M/32 
simulation is 1000M before the peak in 04. Also shown is the 
frequency from the hf = M/32 resolution before the eccen- 
tricity was removed. 

error evident in Fig. 7 is a slowing of the hardening rate 
at low resolutions, causing the final merger to be delayed. 
This delay appears to converge at fourth-order in resolu- 
tion. Viewing the eccentric modulation as a small per- 
turbation on the dynamics of an optimally non-eccentric 
inspiral, we expect that these trends would also provide 
a good approximation for the frequency evolution of a 
merger simulation begun with optimally non-eccentric 
initial data. 


C. Waveform accuracy 

In this section we consider the accuracy of the wave- 
forms generated by our simulations, focusing primarily 
on the accuracy of waveform phasing information in the 
late-inspiral portion of our simulations. Over the course 
of many developmental simulations leading up to these 
results, we found that this early low- frequency part of 
the simulations is the most difficult to accurately sim- 
ulate. This makes sense because the crucial dynamical 
details are in the slow loss of energy and angular mo- 
mentum to the relatively fine, evolving structure of the 
spacetime curvature which ultimately comprises the ra- 
diation. Timescales are also longer for this part of the 
dynamics, requiring high accuracy over a large number 
of computational iterations. 

Fig. 8 compares the 04 waveforms generated by our 
simulations at different resolutions. 

These waveforms have been aligned to coincide at the 
peak in 04 for each of the simulations, as described in 
Ref. [14]. We show the waveform logarithmically because 
the amplitude changes by more than two orders of magni- 
tude through a run, part of what makes these simulations 
challenging. We are especially interested in the phase 



FIG. 8: Waveform at three resolutions. They have been 
shifted in time and phase to agree at the peak of the wave. 
The agreement persists backwards in time, with the growing 
discrepancy in phase converging away with resolution. 


agreement among the different simulations, which is evi- 
dent in the nearly vertical parts of the curves (which ap- 
proach zero crossings) Aligned this way, waveforms gen- 
erated at all three resolutions are nearly exactly super- 
posed after t 250 M. The two higher-resolution simu- 

lations are nearly identical after t ~ — 500M and agree to 
within a small part of a cycle for the full portion of the 
simulation after the initial transient period in the first 
100M of each run. On the other hand, the waveforms 
are easily distinguished by large differences in the overall 
timing of each run, with the lowest resolution simula- 
tion going through a full additional orbit before merger. 
The high-frequency noise evident in the first part of the 
simulations seems to be caused by reflections from our 
refinement boundary interfaces. 

We get a more direct view of phasing information by 
examining the waveforms in polar decomposition. In 
Fig. 9 we show the polarization phase derived from the 
strain rate waveform. 

This time we have aligned the simulations in time from 
the beginning, as should be appropriate for simulations 
of the same initial configuration. Later in the simula^- 
tions, as the frequency increases the phase increases more 
rapidly, and the timing differences among the simulations 
lead to large phase differences. 

We quantitatively compare the phasing results in 
Fig. 10, showing the difference between the two higher- 
resolution runs compared with scaled versions of the dif- 
ference between the two lower resolution runs. We can 
now clearly see that the phase errors (measured this way) 
grow strongly in time. The lower-resolution difference 
has been rescaled two ways, such that they would be ex- 
pected to agree with the higher resolution difference for 
second- and fourth-order convergence, respectively. The 
comparison suggests phasing convergence somewhere be- 
tween second- and fourth-order over much of the run. 
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FIG. 9: Strain rate phase. The high-resolution simulation 
goes through about 14 cycles before merger. The lower reso- 
lution simulations take longer to merge, and go through more 
cycles. 



FIG. 10: Strain rate phase three-point convergence. The 
higher resolution phase difference is shown with and without 
a phase shift to allow comparison of errors at similar dynam- 
ical points in the simulation. After shifting, the phase-error 
appears to be fourth-order convergent for much of the simu- 
lation. 


Since the phase error grows strongly in time, it is inter- 
esting to consider the effect of timing errors in these con- 
vergence comparisons. The first curve in the plot shows 
the high-medium difference shifted in time by 5 7M, rep- 
resenting the timing difference between the two higher 
resolution runs late in the simulation. Viewed this way 
we see time-aligned phase differences taken from similar 
points in the physical evolution of the medium and high 
resolution runs. As a result, the runs appear to converge 
at a faster rate, closer to fourth-order, indicating that 
the timing differences can have a significant impact on 
convergence estimates. 

Above, we have compared our simulations made at dif- 
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FIG. 11: Frequency-based phase comparisons for runs at three 
resolutions. The difference between the lower resolution runs 
is scaled so that it should superpose with the higher-resolution 
difference for fourth-order convergent phasing. The dashed 
curve shows a fit to lj ~ 5 phasing error, which might be ex- 
pected if there are energy conservation violations which are 
constant in time. 


ferent resolutions by comparing the waveforms at equal 
points in time, with time aligned either at the beginning 
of the simulation (t = 0 in the original run) or at the peak 
in -04. Since errors cause the simulations to accelerate 
through their dynamical processes at somewhat different 
rates, such comparisons end up relating moments of quite 
different dynamics, and become less meaningful, depend- 
ing significantly on the reference time according to which 
the phases are compared. This sort of comparison would 
not be applicable at all when there is no clear way to 
physically align predictions in time, such as when com- 
paring numerical simulations with post-Newtonian (PN) 
results. We can avoid assigning a reference time by using 
the gravitational- wave frequency as the reference. 

For the case studied here, a quasi-circular inspiral of 
comparable-mass black holes, it is appropriate to con- 
sider the waveform frequency u = d<j)/dt as a gradually 
increasing monotonic function of time, lj = Though 
the actual waveform frequency of our simulations is not 
monotonic because of small eccentricities in the inspiral 
trajectories, we have shown that we can fit out the ec- 
centric deviations. The residual secular trend u> c (t) pro- 
vides a monotonic frequency which we can apply as an 
independent variable against which to compare various 
simulations. 

We show frequency- based phase comparisons among 
the simulations in Fig. 11, which allows us to compare 
the differences between the two higher resolution simula^- 
tions against the differences between the two lower res- 
olutions. The phase differences are meaningful only up 
to an overall constant; this freedom amounts to applying 
a one-time rotation to a simulation. We have translated 
the phase differences vertically so that they vanish at 
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u = 0.05423, the frequency 1000M before merger in our 
high-resolution simulation. If the simulation errors are 
fourth-order convergent, then the lower-resolution dif- 
ference should be approximately 2.784 times the high- 
resolution difference, as is evident in the figure. 

As is also clear in the figure, phase error accumulates 
most significantly in the low-frequency portion of the sim- 
ulation, Mu; < 0.08. This makes sense generally since the 
simulation spends much more time at lower frequencies. 
Consider, for instance, the effect of a small nonconser- 
vative leakage of energy SE from the simulation. When 
the black holes are well separated, the dynamical devel- 
opment manifested in the sweeping frequency is driven 
by a slow loss of energy and angular momentum. An en- 
ergy leakage SE will change the frequency sweep-rate u 
by SCj ~ (SE)/(dE/duj)j where dE/duj indicates how the 
binding energy changes with frequency. Phase error 5<fi 
can be determined from the error in the sweep-rate by 

5<j> = S J ^-dw = - J ^ Swduj (1) 

Applying the leading-PN-order expressions for dE/ du 
and u(uj) gives a result proportional to lo~ 5 SE. Indeed 
this dependence fits our phase differences rather well, as 
shown in Fig. 11. Note, however, that this does not sin- 
gle out energy nonconservation as the source of error, 
as other errors may produce similar effects. A similar 
leakage of angular momentum would lead to error pro- 
portional to lj~ 4 which also fits reasonably well. 

The clear fourth-order convergence seen in Fig. 11 sug- 
gests that we can apply Richardson extrapolation to es- 
timate the difference of our high-resolution (M/ 32) run 
from the infinite resolution limit. This provides a phase 
error estimate for this run of 0.93 times the difference 
between the phases from the 3M/80 and M/32 resolu- 
tion runs. Note that for the last 1000M of our M/32 
simulation, we estimate roughly two radians of phase er- 
ror, which is about a third of a gravitational wave cy- 
cle. A benchmark for accumulated waveform phasing 
errors is one-half of a cycle, because phase error exceed- 
ing this amount would lead to destructive interference in 
matched-filtering applications. For our high-resolution 
simulation, we estimate less that one-half cycle of grav- 
itational wave phase error over the full simulated wave- 
form, excluding the meaningless transients in the first 
100M. In [22], we estimated that these phasing errors 
are smaller than the implicit phasing difference between 
the 3PN and 3.5PN expansions of u(u) after t ~ — 300M 
(c u c M ~ 0.08). 

Frequency-based phase comparisons, such as we have 
presented here, are better suited than time-based phase 
differences, which depend strongly on where the wave- 
forms are chosen to be aligned in time. The relationship 
between the two can be understood by considering a one- 
parameter family of waveform results, with a parameter 
A representing model dependence, in this case the numer- 
ical grid-spacing. The waveforms would provide phase as 


a function of time from which we can derive fre- 
quency which is monotonic for small eccentricity. 

Inverting to obtain t\(u) one can derive the frequency- 
based phasing 4>\{u) = <t>\{t\{w)). Now considering vari- 
ations S = d/dX near A = 0, one finds the relationship 
between frequency- and time-based phase comparisons 

54>{uj) = 6(f>(t(u)) + u>St(u>). (2) 

This sheds some light on the often confusing issue of time 
alignment in the time-based comparisons shown earlier. 
Specifically, for a waveform which sweeps significantly 
through frequency, time- and frequency-based phase dif- 
ferences will be most similar when the time-based phase 
differences are aligned so that St vanishes where ui is 
largest. In our case, the net frequency-based phase dif- 
ferences in Fig. 11 are closer to net phase differences with 
time aligned at the end of the waveform, as in Fig. 8, than 
when time is aligned at the lower-frequency beginning, as 
in Figs. 9 - 10. 

IV. APPLICATIONS TO GRAVITATIONAL 
WAVE OBSERVATIONS OF BBHS 

In the first part of this paper, we presented a state-of- 
the-art calculation of the late inspiral and merger of equal 
mass, nonspinning BBHs, starting ~ 7 orbits before the 
peak radiation amplitude. In Ref. [22] we showed that 
there is a significant region over which the waveform from 
the M/32 numerical simulation presented above agrees 
with waveforms derived from PN calculations. In this 
section we will utilize the best of both treatments by 
joining the final segment of the BBH evolution with the 
PN approximation for the preceding inspiral, as shown 
in Fig. 12. We will show below that the phasing infor- 
mation for this waveform has an estimated accuracy of 
better than one-half of a gravitational wave cycle, making 
the waveform applicable in a variety of matched-filtering- 
based gravitational wave data analysis applications. We 
use this waveform to calculate the relative detectability 
of the inspiral and the merger-ringdown, as well as the 
detectability of the entire waveform, for iLIGO, adLIGO, 
and LISA. We show examples of characteristic signal 
strains for both classes of detectors. We also compute 
SNRs for astrophysically interesting BBH masses, high- 
lighting the importance of the late inspiral and merger 
parts of the signal. Modifications of our results that may 
arise from BBHs with spins and unequal masses are con- 
sidered. 


A. Waveform matching and SNR calculation 

The PN waveform used for all the analyses in this sec- 
tion is of order 3.5PN in phase (as derived in [37, 38]) 
and 2.5PN beyond the quadrupole approximation in am- 
plitude (as derived in [39]). It should be noted that the 
PN approximation that we use for the phase is actually 



10 



m 


FIG. 12: Numerical simulation results for gravitational wave strain compared with post-Newtonian estimates. The waveform 
shown is from the high resolution numerical simulation presented in Sec. II overlaid here with a PN waveform with 3.5PN order 
phasing and 2.5PN order amplitude accuracy. The combined waveform, joined at t = —32 8M (circle) is applied in Sec. IV to 
calculate signal- to-noise ratios for iLIGO, adLIGO, and LISA. 


an expansion of the chirp rate u in terms of the frequency 
u, which we then numerically integrate. Direct numeri- 
cal integration of the 3.5PN expansion of the chirp rate, 
V 3 . 5 PN, for example in the integrand du(u /u> 3 . 5 pat), does 
not strictly respect the PN approximation of the phase, 
as the latter would require additional 3.5PN expansion of 
the integrand itself. However, the phase obtained in this 
manner will have the same convergence properties as the 
original u^^pn expansion, which is arguably a more fun- 
damental PN quantity because of its close relationship to 
the rate of energy loss, E, from which it is derived. Ad- 
ditional expansion of the phase appears to compromise 
its accuracy. Using shorter runs it had been observed in 
[20] that the phase obtained from numerical integration 
of the PN expansion of the chirp rate seemed to agree 
better with numerical simulations than did the analytic 
expressions; more recently it was demonstrated in [22], 
using runs extending well into the late inspiral and with 
the effects of eccentricity mitigated, that the numerically 
integrated phase appears to be converging to the numer- 
ical result at frequencies where the full 3.5PN expansion 
of the phase is clearly invalid. We therefore use the nu- 
merically integrated 3.5PN chirp rate for the phasing of 
the PN portion of our waveform. 

Fig. 12 shows our numerical waveform from Sec II over- 
laid with the PN waveform that was just described. To 
generate a complete, mass-scalable waveform, we match 
the frequency of the numerical simulation to the PN pre- 
diction, adjusting the phases to also be equal at that 
point as shown in Fig. 12, and connecting the two halves 
to make a single waveform. This is done by shifting the 
PN waveform until the frequency equals the numerically- 
predicted frequency at a time in the simulation where 
the accuracy of the numerical data first surpasses the 
accuracy of the PN approximation, as estimated in [22]. 
Specifically, [22] predicts this point of equivalent accuracy 


to occur at Mu ~ 0.08, which corresponds to t = — 328M 
(shown by the circle in Fig. 12). It is worth noting that 
there was no need to adjust the PN amplitude for con- 
tinuity. The amplitude agreement with the numerical 
simulation is so good, and hence the resulting amplitude 
is so nearly continuous, that the small discontinuity fails 
to produce any discernible artifacts in the Fourier trans- 
form h(f) of the resulting waveform. 

Having generated a waveform, it is informative to esti- 
mate the waveform’s phasing accuracy over the course of 
the BBH evolution. Note in Fig. 11 that for the portion 
of our M/32 simulation that is used in the waveform, we 
estimate <0.5 radians of phase error. If we take the dif- 
ference between 3 and 3.5 PN terms to be an estimate of 
the phase error as in [22], we can assess the error for the 
PN portion of the waveform. It was shown in [40] that 
the analytic PN phase expression accumulates very little 
error, on the order of 0.1 radians, until Mu ~ 1 x 10 -4 . 
Beginning our numerical phase integration at this point 
and evaluating up to Mu = 0.08 yields a phase error of 
~ 1.8 radians, such that the total accumulated phase er- 
ror over the entire waveform is £ 2.5 radians. As stated 
previously, an accumulated waveform phasing error of 
less than one-half of a cycle is the threshold below which 
wave-matching comparisons may be used for matched- 
filtering applications. We therefore have a waveform with 
sufficient accuracy throughout its entire evolution to be 
useful in gravitational wave data analysis. 

The calculated waveform on which we have just de- 
scribed is actually the total strain on the equatorial 
plane, where h x vanishes and therefore h+ provides the 
only contribution. To get the optimal strain amplitude, 
we multiply the result by 2\/2, which is the ratio of peak 
gravitational wave amplitude (which occurs on the equa- 
torial axis) to the amplitude of h+ alone. We then sim- 
ply divide by y/E in order to convert to an orientation- 
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averaged waveform for our subsequent analyses. This 
factor can be understood by observing that orientation- 
averaging is fully equivalent to averaging over all sky po- 
sitions of the detector from the perspective of the BBH, 
and such sky- averaging results in a factor 1 / y/5 in sensi- 
tivity [1]. 

The SNR is calculated assuming matched-filtering is 
performed on the data, and that the waveforms are per- 
fect copies of the embedded signal. In this case, the sky- 
and waveform-polarization-averaged SNR is given by 

< (SNR) 2 >= J d(\nf) (3) 

where h char (f ) = 2/|h(/)| is the characteristic signal 

strain and h n (f) = V5h rms (f) = \/5fS n (f) is the rms 
of the detector noise^ fluctuations multiplied by \/5 for 
sky-averaging, with h(f) and S n (f) being the Fourier 
transform of the signal strain and the power spectral den- 
sity of the detector noise, respectively [1]. 

The waveform scales with luminosity distance Dl and 
total mass M as /i c ^ar oc (1 + z)M / Dl> while the time 
axis for an observed wave, after redshifting, scales as 
t oc (1 -f z)M, so that the waveform shown in Fig. 12 
is applicable over all total masses and redshifts. When 
needed we relate luminosity distance to redshift z using 
cosmological parameters consistent with the most recent 
WMAP results (fbv = 1 - fijw = 0.72, h = 0.73) [41] and 
the relation 


Pi(2) = ii±£)£ r dz ' 

Ho Jo 1 + z ') 3 -j- Da 


(4) 


For cases where an impractically long time series would 
be needed to cover the band with an adequate sampling 
rate, the waveform is extended in Fourier space to still 
lower frequencies (and consequently back further in time) 
using the quadrupole formula, 


^quad(/) — 
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(5) 


The PN portion of the waveform continues slightly past 
where its Fourier transform deviates from the quadrupole 
expression for h(f) by ~ 2%, at which point (5) is used 
to extend h(f) back as far as necessary. The PN seg- 
ment is truncated at a higher frequency than the low- 
est frequency component of its Fourier transform in or- 
der to eliminate edge effects. Finally, the quasi-normal 
ringdown at the end of the numerical simulation is ex- 
tended by fitting a damping coefficient and fundamental 
frequency to the data in order to mitigate edge effects at 
the high-frequency end of the Fourier transform. 


B. Observing Stellar BBHs and IMBBHs with 
iLIGO and adLIGO 


Ground-based interferometers are sensitive to rela- 
tively high frequency gravitational waves from coalesc- 
ing stellar mass (M < 10 2 M G ) and intermediate mass 
(IM) (10 2 M g < M < 10 3 M 0 ) BBHs. In this section, we 
apply our combined waveform to consider the response 
of iLIGO and adLIGO to BBH coalescence, illustrat- 
ing the importance of numerical simulation results for 
ground-based detectors. LIGO has facilities in Hanford, 
WA and Livingston, LA; each facility has an interfer- 
ometer consisting of two 4 km long arms, and Hanford 
also has a 2 km detector. Both the iLIGO and adLIGO 
detectors are designed to detect high frequency grav- 
itational waves, sith iLIGO sensitive in the frequency 
range 40Hz < / < 8000Hz and adLIGO in the range 
14Hz < / < 10 3 Hz. Initial LIGO is currently operating 
at or near the initial design sensitivity in a year-long sci- 
entific data^taking run. Advanced LIGO is a planned 
upgrade that will increase the detector sensitivity by 
roughly an order of magnitude across the band. In addi- 
tion, adLIGO can be tuned to optimize its sensitivity for 
different sources. 

For our analysis of iLIGO, we used the design sensi- 
tivity to characterize the detector noise [42]. This sensi- 
tivity assumes that the noise is seismically limited below 
40 Hz, thermally limited between 40 and 150 Hz, and 
shot-limited above 150 Hz. For adLIGO, unlike iLIGO, 
we had a choice of tuning configurations. We used the 
wide-band tuning typically associated with burst sources 
because of its dramatically superior sensitivity at higher 
frequencies, where the merger portion from many sources 
is predicted to occur [43]. This yielded an improved SNR 
for most masses compared to tunings that were optimized 
for only the early inspiral portion of the coalescence. 

In Fig. 13, we show h c har for several sources plotted 
relative to the h rms sensitivity curves for iLIGO (dashed 
line) and adLIGO (dash-dotted line). We plot these val- 
ues because the height of h c har above h rms is an indicator 
of the SNR, as can be seen by inspecting Eq. 3. 

By rescaling we can calculate the SNR as a function 
of redshifted mass, and particular luminosity distance 
Dl. In Fig. 14 we plot the SNR achievable by iLIGO for 
sources at a luminosity distance Dl = 100 Mpc as a func- 
tion of redshifted mass (1 + z)M. Here, the dashed line 
shows the SNR from the early inspiral in the time range 
— oo <t< —1000 M, which is roughly up to the start of 
our run. The dotted line shows the SNR for the late inspi- 
ral, — 1000M < t < —50 M, where t = —50 M is approx- 
imately the time at which the merger burst begins. The 
thin solid line gives the SNR for —50 M <t< oo, and en- 
compasses the merger-ringdown part of the signal. The 
thick solid line shows the SNR from the entire waveform. 
Note that the addition of the merger-ringdown waveforms 
increases the SNR and extends the detectable mass range 
significantly. The merger-ringdown portion t > —50 M 
dominates for all equal-mass nonspining merger observa- 
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FIG. 13: The iLIGO (dashed) and adLIGO (dash-dotted) rms 
noise amplitudes h n with the characteristic amplitudes h c har 
of 6 example sources (solid). The locations on each h c h ar 
corresponding to the peak 04 amplitude (circle) and 1 second 
before the peak in the observer’s frame (filled circle), as well 
as t = —50 M (square) and t = — 1000M (diamond) in the 
source’s frame, are as marked. 



FIG. 14: SNR for sources at luminosity distance Dl — 100 
Mpc plotted vs. redshifted mass for iLIGO. The contributions 
from -oo < t < — 1000M (dashed), -1000 < t < —50 M 
(dotted), and —50 M < t < oo (thinner solid), as well as the 
SNR from the entire waveform (thicker solid) are shown. 


tions detectable with SNR larger than 10 at 100 Mpc. 

This type of plot was first made in Ref. [1], and it 
is useful to compare our results with theirs. Our SNR 
calculations are based on a full waveform for the case of 
equal-mass, nonspinning black holes. The work in Ref. 
[1] was done before merger waveforms were calculated 
and thus is based on estimates for the merger-ringdown 



FIG. 15: SNR for sources at luminosity distance Dl = 1 Gpc 
plotted vs. redshifted mass for adLIGO. The contributions 
from — oo < t < — 1000M (dashed), —1000 < t < -50 M 
(dotted), and —50 M < t < oo (thinner solid), as well as the 
SNR from the entire waveform (thicker solid) are shown. 


regime. For example, they estimated a merger radiation 
efficiency of ~ 10%, which is higher than our results but 
may well obtain for mergers with spin. Comparing their 
Fig. 4 for the SNR for iLIGO with our Fig. 14 we note 
that their curve for the inspiral includes the radiation up 
to the merger and so should be compared to the combi- 
nation (in quadrature) of our dashed and dotted curves. 
Our result for the merger SNR is somewhat smaller than 
theirs, due to the smaller amount of radiation emitted in 
our mergers. More recently, an analysis of SNR for iLIGO 
using numerical relativity waveforms for the merger and 
PN waveforms for the inspiral was made in Ref. [20]; our 
results in Fig. 14 are similar to what they report in their 
Fig. 22. 

Fig. 15 shows the SNR for sources at Dl = 1 Gpc for 
adLIGO. Comparing with Fig. 14, we see that adLIGO 
will have a significantly higher sensitivity to BBHs over 
iLIGO. This point is reinforced in Fig. 16, which shows 
contours of SNR for adLIGO as functions of redshift 
z and total mass M. We find that for M ~ 200 M©, 
adLIGO should be able to achieve an SNR greater than 
10 out to nearly z = 1 for equal mass nonspinning bina- 
ries. From Fig. 14 it is evident that these high SNR’s de- 
pend strongly on the merger-ringdown part of the wave- 
form t > —50 M. 

It is important to note that astrophysical BBHs are 
likely to have mass ratios different from unity, and that 
this will reduce the SNRs computed here for the equal 
mass case. For stellar BBHs, current work [44] shows 
mass ratios mi/rri 2 ~ 0.8 [45] however, this is based on 
small number statistics and may change when larger sim- 
ulations are completed. The rates for such mergers may 
be low, ~ 2 yr~ l for adLIGO, depending on the evolu- 
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FIG. 16: SNR contour plot with mass and redshift depen- 
dence for adLIGO. 


tion of the original binary through the common enve- 
lope phase [46]. For IMBBHs, mass ratios in the range 
0.1 < mi/rri 2 < 1 are expected to be the most relevant, 
with potential rates of ~ 10s per year [47]. We can ap- 
ply the mass scalings from Ref. [1] to show the effect 
of mass ratios on the computed SNRs; specifically, SNR 
~ r? 1 / 2 for the inspiral, and SNR ~ 77 for the merger and 
ringdown, where 77 = fi/M and (i = 77 a 1 7722 /M is the re- 
duced mass. Astrophysical BBHs are also expected to 
be spinning and this can potentially raise the SNR, for 
example if there is a spin- induced hangup that generates 
more gravitational wave cycles in the merger [18]. 


C. Observing MBBHs with LISA 

LISA, a proposed space-based interferometer consist- 
ing of three 5 million km long arms, will be sensi- 
tive to low-frequency gravitational waves from coalescing 
MBBHs in the band 3 x 10 -5 Hz < / < 1Hz. The coalesc- 
ing massive binary black holes (MBBHs) which radiate 
in this band will have masses M > 10 4 M©. 

Fig. 17 shows h C har for several MBBHs plotted rela- 
tive to the LISA sensitivity curve. We used the “stan- 
dard” LISA sensitivity curve [48, 49] for frequencies 
above 1 x 10 -4 Hz, with shot and pointing noise con- 
tributions totaling 20pm/\/Hz of laser phase noise. For 
3 x 10 _5 Hz </< lx 10~ 4 Hz, we employed a more con- 
servative estimate of the acceleration noise than the one 
given in [48] , instead assuming a steeper amplitude spec- 
tral density which falls off as /~ 3 constrained to match 
the standard sensitivity curve at 1 x 10~ 4 Hz [50]. Below 
3 x 10 _5 Hz, we assume the detector has no sensitivity, 
which is a reflection of the uncertainty of the sensitivity 
at such low frequencies and our desire to make conserva- 
tive estimates. The sensitivity model assumes that there 



FIG. 17: LISA rms noise amplitude firms from the detector 
only (dashed) and from the detector combined with the antic- 
ipated white dwarf binary confusion (dash-dotted) [51] with 
the characteristic amplitudes h C har of three example sources 
(solid). The locations on each h c har curve corresponding to 
the peak amplitude (circle), 1 hour before the peak (filled 
circle), 1 day before the peak (circle with inscribed cross), 
and 1 month before the peak (circle with inscribed square) 
in the observer’s frame, as well a s t = —50 M (square) and 
t = — 1000M (diamond) in the source’s frame, are as marked. 
When necessary, the quadrupole approximation is used to ex- 
tend hchar backward in time 3 years before the peak 7/4 am- 
plitude in the detector’s frame (dotted). 


are no correlated noise sources, and its power character- 
ization corresponds to a round trip through one arm of 
the interferometer. 

MBBH sources can remain in-band for LISA over a 
very broad frequency range. Therefore, unlike the case 
of iLIGO and adLIGO, LISA sources nearly always re- 
quire the use of the quadrupole approximation procedure 
mentioned above to extend h C har to sufficiently low fre- 
quencies. Also, since more massive BBHs chirp more 
slowly, a MBBH could potentially be in LISA’s sensitive 
band for much longer than the mission’s lifetime. To 
prevent unrealistic SNR values due to this excessive in- 
tegration time, the quadrupole formula is only used to 
extend h c har to a low enough frequency such that the to- 
tal h c har used in our calculations corresponds to 3 years 
of data in the detector’s frame, which is a conservative 
estimate of the expected mission lifetime. 

The SNR for LISA is shown as a function of redshifted 
mass, normalized for Dl = lOGpc in Fig. 18. The bump 
in the curves is caused by the binary confusion noise. 
Again we see the enhancement of SNR from the merger- 
ringdown part (thin solid line) of the waveforms, and 
confirm the strikingly large values of SNR obtainable by 
LISA for these sources seen in [1] and [20]. For systems 
with redshifted mass (1 -F z)M < 3 x 10 4 , the early inspi- 
ral t < — 1000M portion of the waveform dominates. The 
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FIG. 18: SNR for sources at luminosity distance Dl — 10 
Gpc plotted vs. redshifted mass for LISA. The contributions 
from the early inspiral — oo < t < — 1000M (dashed), late 
inspiral —1000 < t < —50 M (dotted), and merger-ringdown 
—50 M < t < oo (thinner solid), as well as the SNR from the 
entire waveform (thicker solid) are shown. 


highest SNRs for equal-mass nonspinning mergers are ob- 
tained for systems with(l + z)M > 10 6 , again dominated 
by the merger-ringdown portion of the waveform. 

Contours of SNR for LISA are shown in Fig. 19 and 
demonstrate that LISA can observe MBBHs through- 
out the observable universe at large SNRs. We find 
it encouraging that, in addition to the large SNR val- 
ues predicted for LISA overall, some of the largest val- 
ues out to the largest redshifts occur in the mass range 
10 5 Mq < M < 10 7 M 0 where models of BBH popu- 
lations predict that the binaries can coalesce within a 
Hubble time [52] and that the event rates for LISA are 
several per year [53]. As discussed above, the effects of 
unequal masses will tend to decrease these SNR values, 
while spins may increase or decrease them. 

Even for non-optimal configurations, the presence of a 
MBBH coalescence in the LISA data stream can domi- 
nate all the anticipated noise sources. Fig. 20 shows a 
simulation of LISA’s response to the merger of equal- 
mass nonspinning black holes with total mass M = 
10 5 M© located at redshift z = 15, and oriented so that 
LISA lies in the system’s equatorial plane, where the ra- 
diation is weakest. The SNR for a signal from such a 
source will be ~ 200 averaged over sky positions and po- 
larizations (see Fig. 19). 


V. SUMMARY AND DISCUSSION 

Coalescing BBHs are expected to be the strongest 
sources for both ground-based interferometric as well 
as the space-based LISA. In particular, the strong-field 



m 


FIG. 19: SNR contour plot with mass and redshift depen- 
dence for LISA. Note that MBBHs with masses M > 10 7 Mq 
may not coalesce within a Hubble time [52]. 



FIG. 20: Simulated LISA data stream showing LISA’s re- 
sponse to a system of two equal mass black holes (M = 
10 5 M©) located at redshift z=15 observed on the system’s 
equatorial plane. The quantity plotted is an unequal arm 
Michelson interferometer observable “X” [54]. The LISA re- 
sponse and instrumental noise are realized using the LISA 
Simulator [55, 56], and colored noise was added to represent 
the unresolvable galactic binary foreground with the spectrum 
used in Ref. [51]. The inset shows the signal over a longer du- 
ration where low frequency noise is evident. 


merger portion of the gravitational wave signal produces 
an intense burst of radiation and has the highest luminos- 
ity, emitting more energy per second than the combined 
starlight emitted in the observable universe. 

Recent breakthroughs in numerical relativity have 
opened a new era in understanding the late stages of bi- 
nary black hole coalescence. We now have a good under- 
standing of the merger-ringdown signal, starting ~ 50 M 
before the peak radiation amplitude, for equal mass non- 


0, (Gpc) 
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spinning BBHs. The late inspiral evolution, that is, more 
than a few orbits prior to ringdown, is more challenging. 
Such simulations require better numerical stability and 
more computational resources, as well as higher accuracy 
to control the accumulated phase error. 

In this paper, we presented new simulations of 
equal mass nonspinning BBHs starting in the late 
inspiral regime and covering approximately an addi- 
tional factor of three in frequency before the merger- 
ringdown. We carried out runs at three resolutions, 
hf = 3M/64,3M/80, and M/32. Our runs start with 
relatively low eccentricity and show good convergence 
and conservation properties. We have demonstrated the 
stability and accuracy of our simulation over the course 
of this unprecedented seven orbits. We also showed the 
value of using frequency (rather than time) to set a refer- 
ence for the purpose of comparing results between runs as 
well as with the PN approximation. In recasting phase 
vs. frequency we have found particularly good agree- 
ment, not only between the runs but also with PN pre- 
dictions. 

We have also matched the resulting gravitational wave- 
forms to PN calculations covering the earlier parts of the 
inspiral. The resulting full waveform has less than half a 
cycle of accumulated phase error over its entire frequency 
band. Using this waveform, we calculated the SNRs for 
iLIGO, adLIGO, and LISA. Our results confirm the im- 
portance of the merger-ringdown signal, which yields the 
highest values of SNR for the majority of equal mass 
signals [1, 20]. We also show the SNR for the late in- 
spiral regime, which numerical simulations are now be- 
ginning to simulate. The late inspiral dominates the 
SNR for iLIGO and adLIGO for the lower mass (< a 
few xlOM©) stellar BBHs, and the SNR for LISA gen- 
erated by MBBHs with M ~ 10 1 2 3 4 5 6 7 8 M 0 . Contour plots of 
SNR as a function of z and M show that adLIGO can 
achieve SNR > 10 for IMBBHs out to nearly z ~ 1, and 
that LISA can observe MBBHs at SNR > 100 out to the 
earliest epochs of structure formation at z > 15. 

Our work has focused on equal mass nonspinning 


BBHs. Astrophysically, BBHs are expected to have un- 
equal masses and spins. In general, the effects of un- 
equal masses will tend to decrease the resulting SNRs, 
while spins can increase them. Calculations of merger- 
ringdown waveforms for several mass ratios, and for some 
spins are currently available; we expect this to be a sig- 
nificant area of focus in the foreseeable future, both ex- 
panding the range of parameters studied and extending 
the duration of the resulting simulations. 

As the computational technologies mature, simulations 
of the merger can be used in conjunction with gravita- 
tional wave observations to probe gravity in the arena 
of strong fields. In particular, if the binary masses and 
spins can be obtained from observations of the inspiral 
(as demonstrated, e.p., in Ref. [57] for LISA), numerical 
relativity can be used to calculate the merger waveform. 
This will allow a comparison between the predictions of 
general relativity - or indeed, any other theory of gravity 
used in a numerical simulation - with observations in the 
regime of very strong gravity. 
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